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Abstract 



For ferrofluids that are described by a system of hard spheres interacting 
via dipolar forces we evaluate the magnetization as a function of the internal 
magnetic field with a Born-Mayer technique and an expansion in the dipolar 
coupling strength. Two different approximations are presented for the mag- 
netization considering different contributions to a series expansion in terms 

of the volume fraction of the particles and the dipolar coupling strength. 
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I. INTRODUCTION 

Ferrofluids (l| are suspensions of ferromagnetic particles of about 10 nm diameter in a 
carrier fluid. The particles are stabilized against aggregation by coating with polymers or by 
electrostatic repulsion of charges brought on their surface. On macroscopic scales, ferrofluids 
can be described as liquids with intrinsic superparamagnetic properties. 

In this paper, we are concerned with the equilibrium magnetization M as a function of 
the internal magnetic field H for given temperature and particle concentration. Sufficiently 
low concentrated ferrofluids behave like a paramagnetic gas. Therein the interaction be- 
tween the particles can be neglected and the equilibrium magnetization can be described 
properly by the Langevin function. The magnetic properties are then necessarily weak. To 
produce ferrofluids with strong magnetic properties one has to have either a higher particle 
concentration or one has to use ferromagnetic material with a large bulk magnetization, 
e. g., cobalt instead of magnetite. In both cases the magnetization is strongly influenced by 
dipole-dipole and other interactions between the particles. 

Several models of dipolar interacting systems have been studied in the literature. Nu- 
merical investigations were based on density functional approaches 0-|J and Monte Carlo 
simulations . The models differ in the treatment of the short range interactions, which 

were described by hard sphere P,P,p|,p!3|-p!5[| , other hard core potentials P,|l2||, soft sphere 
|7|JlO[|, or Lenard- Jones potentials |F2|j3|j6|,f?|,p]1 . These investigations were mainly undertaken 



to reveal the phase transition properties. These properties are substantially different for 
different short range interactions. Thus, for example the question whether a system of par- 
ticles interacting via long range dipolar forces shows without any dispersive energy, e. g., 
from attractive van der Waals energy a "liquid-vapor" phase coexistence of a dense and a 
less dense phase is currently being discussed [[f|, p!2|p!4| -[T8[| . 

Analytical models focus mainly on the equilibrium magnetization in the gas phase (were 
the term "gas" refers, as far as ferrofluids are addressed, to the magnetic particle subsystem 



within the liquid carrier). Such models are the Onsager model RI5f|, the Weiss model [GO], the 
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mean spherical approximation [^T] , and an approach by Buyevich and Ivanov J22| (called high 



temperature approximation in |23j|). These models were tested experimentally for ferrofluids 
|23^f26|. Especially the mean spherical model and the high temperature approximation 
showed good results P^| . 



Our approach assumes the magnetic particles in the ferrofluid to be hard spheres with 
a common diameter D and dipolar moment m. We use the technique of the Born-Mayer 
expansion |27| together with an expansion in the strength of the dipolar coupling to get 
analytical approximations. They are obtained via series expansions of the free energy in 
terms of two parameters: (i) the volume fraction of the hard core particles and (ii) a 
dimensionless dipolar coupling constant e, given by the ratio between a typical dipolar 
energy for particles in hard core contact and the thermal energy kT. Our result for the 



magnetization goes beyond the high temperature approximation [2^] and reduces to it in 
linear order in and e. 

Dipolar forces fall off as r -3 and are thus of long range nature. This long range character 
requires great care when invoking the thermodynamic limit p8|-|30|j. To circumvent the 
problem we model the dipolar fields that are generated by distant particles by a magnetic 
continuum field (similar to the treatment in the Weiss model) while incorporating the near 
field contributions explicitly in a statistical mechanical description. The magnetization M 
is then derived as a function of the internal magnetic field H. The so obtained relation 
M(H) is independent of the probe geometry. Once M(H) is known, the magnetization for 
a given geometry can in principle be derived by solving the macroscopic Maxwell equations. 
This may practically still be a difficult task, at least as long as the external field is small or 
absent. In this case it is known that the magnetization will show for general shaped probes 
a nontrivial spatial variation at high enough densities . 

Since our method yields an expression for the free energy of the model system we can 
in principle calculate also other thermodynamic quantities and in particular address the 
questions of phase transition, e. g., between gas and liquid or between ferromagnetic and 
non-ferromagnetic phases. We have not addressed the question of a gas-liquid transition of 
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the magnetic particles suspended in the ferrofluid since it is believed that short range van der 
Waals-like attractions would have to be incorporated to model real ferrofluids appropriately 



in this regard [Hj. However, the question whether a strong dipolar coupling induces in 
zero external field a spontaneous magnetization that is currently debated in the literature 
PJ3|]^,P, |31| , |3^ 1 is briefly touched upon in Sec. |V1 A| of this paper. 

The paper is organized as follows: In Sec. |T| we discuss the connections between the 
various fields that are of relevance in a ferrofluid. We present the model to treat the long 



range dipolar forces. In Sec. [HJ we present the expansion method to get analytical solutions 



in terms of the two small parameters e and 0. In Sec. [TV] we calculate an expression for the 
magnetization that contains only linear terms in <p but, at least in principle, arbitrary high 
orders in e. In Sec. [V| a different expression is derived containing also quadratic terms in 



but also only up to second order terms in e. In Sec. VI we discuss our findings and investigate 



the applicability of the results in the 0-e plane. Sec. VII contains a short conclusion. 



II. MAGNETIC FIELDS AND MAGNETIZATION 

We are interested in the effect of dipolar interactions of the magnetic particles in a 
ferrofluid on the equilibrium magnetization of the ferrofluid. To that end we consider the 
ferrofluid as an ensemble of identical spherical particles of diameter D, each carrying a 
magnetic moment of magnitude m. These particles interact with each other via magnetic 
dipole-dipole interactions and a hard core repulsion with hard core diameter D. We assume 
D > D mag where D mag is the diameter of the magnetic core of the particles thus allowing 
for a surfactant surface layer that provides a steric repulsion. 

The particles can lower their potential energy by orienting their magnetic moments par- 
allel to a local magnetic field. However any interaction of the particles with the fluid medium 
they are suspended in is ignored. The latter is taken to be magnetically inert. 
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A. Different magnetic fields 

Before we outline in Sec. [II B| how we determine in principle the magnetization of the 
ferrofluid we should like to review briefly the different magnetic fields that one has to dis- 
tinguish and that are of importance in a system with dipolar interactions. The first field is 
the external magnetic field H e that is applied outside of the probe. If dipolar interactions 
can be neglected, H e is also the field at the position of the particles - at least as long as 
the carrier fluid can be treated a magnetic vacuum which we will assume throughout the 
paper. In the presence of dipolar interactions additional fields have to be considered. One 
of them is the internal field H, that is the macroscopic field inside the probe. By assuming 
that the equation of state M(H) is known, H can be calculated by the common methods 
of continuum magnetostatics. But the macroscopic field H differs in general from the field 
H; oca ; that the magnetic particles feel. 

So far one has employed in the ferrofluid literature two models to calculate ttiocai from 
H that are similar in spirit, namely the Weiss model |20] and the Onsager model |0J. Both 
introduce a virtual hollow sphere inside a magnetic continuum such that the sphere contains 
a single magnetic particle in its center. The Weiss model assumes the magnetization M and 
the internal field H to be constant everywhere in the magnetic continuum surrounding the 
sphere. Then the field inside the sphere is given by 

M 

H s = H+__. (2.1) 
This is the field that the single magnetic particle feels within the Weiss model, i. e. Hi oca i = 

The Onsager model, on the other hand, is restricted to linearly responding fluids and 
calculates the field inside the sphere on the assumption that it is really hollow and that 
therefore H and M differ near the sphere from its bulk values. In that case the field within 
the sphere is 

Hiocai = n X , H . (2.2) 
2y + 3 
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X is the susceptibility. In both models the magnetization is calculated as the magnetization 
of a system of noninteracting dipoles in the magnetic field Hi OC ai, i- e., 

M = M sat C ( —Hi oca i ) . (2.3) 



Here 



A" m /n , s 

M sat = —— 2.4 
V Ho 

is the saturation magnetization of the fluid, £ the Langevin function, m the magnetic 
moment of the particles, and N/V their number density. In the Onsager model the Langevin 
function is consistently used only in linear order. Letting M = %H. on the left hand side of 
( P- 3D and using ( |2.2| ) allows to calculate x- 

In the Weiss model the selfconsistent solution M(H) is determined using ( |2.3p and ( |2.1| ). 
The Onsager and Weiss models differ in the treatment of the back reaction of the particle 
inside the sphere on the magnetic continuum near the sphere's boundary. 



B. Decomposition of fields 

Since the magnetic continuum is a macroscopic concept one should be careful when 
using it on the mesoscopic length scales of the interparticle distances and particle diameters. 
A first -principle statistical mechanical calculation of the magnetization would start with 
expressing the energy of the system in terms of the statistical variables of the constituents. 
In this context the local magnetic field H/ oca / that a magnetic moment, say, at position 
Xj feels is of importance. It is composed of two different magnetic fields, the external 
field H e and the dipolar contribution Hdipoie from the other — 1 particles at positions 
Xj, possessing a magnetic moment nx,. Thus within the first-principles approach one has 
Hiocai = H e + H dipo i e , where 

■ 47T/i r i:; - 

Here r y x; Xj , r^j | r ^ | , and r ^ r ^ jv^j . 
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The long range character of the dipolar forces requires special care when invoking the 
thermodynamic limit V — > oo and iV — > oo (N/V = const.) — for a critical discussion see, 
e.g., Ref. [^,^]. The reason is that the dipolar contribution ( |2.5|) will in general depend on 
the geometry of the ferroffuid probe and the location of Xi within it. Thus the equilibrium 
magnetization of a probe in an external field will in general depend on the geometry of the 
latter and furthermore it will be spatially varying. We therefore use here an approach similar 
to the one that has been used successfully in solid state theory |$3[ to determine, e. g., the 



crystal field splitting caused by local fields. It properly accounts for the contributions from 
microscopic and macroscopic scales. 

Consider some magnetic particle i in a ferroffuid probe in thermodynamic equilibrium. 
The particles beyond some distance R s from Xj can be considered as independent from 
particle % if R s is larger than the correlation length induced by the dipolar interactions. 
Furthermore, if R s is large enough, their contribution to the local field at Xj can be approxi- 
mated by a contribution from a magnetic continuum with equilibrium magnetization M and 
macroscopic field H. We assume that the distance R s is still small compared to the length 
scale on which the macroscopic fields M and H vary. Thus we introduce a virtual sphere 
of radius R s around particle % (dark particle in the center of Fig. |TJ) to separate the dipolar 
field into a "far" and a "near" contribution 

■ I / \ _ \p 3f,;j(mj ■ ijj) — nij 

"-dipole.for V- x -v / v a % > \ / 

II / \ _ v-> 3r^(mj ■ ijj) — m ; 

Then 

H; oca ; H e -\- Hdipole,far ~Hdipole,near • (^"^) 

If the sphere would be empty, H e + Hdipoiejar = H s would be the local field inside the sphere. 

A key point of our treatment is to express the far field ttdipoiejar within the continuum 
approximation. Using this approach the field in the empty sphere is given by H s = H + M/3 



February 1, 2008 8 



fl2.ip . Note that this approximation is not valid near the sphere's boundary. But at the center 
of the sphere H; oca / consists of the continuum contribution H s from the far region and the 
contribution for the dipoles within the sphere: 

M 

Hlocal H s -|- Hdipole,near H -|- ~\~ Hdipole,near ■ (2'9) 

This result agrees with Eq. (27.26) of Ref. |33 where it has been derived with slightly 



different arguments for electric dipoles. Due to the long range character of the dipolar 
forces Hdipde will be in general geometry dependent and spatially varying. H s , respectively 
H and M will then also show these features as mentioned above. 

If the dipolar coupling between the particles is so weak that even the dipolar fields of 
the nearest neighbours of particle i can be described by a continuum field, i. e., if R s can 
be chosen as being smaller than the mean distance between the particles, we can drop the 
contribution Hdipoie,near altogether and arrive at the Weiss model, where a single particle is 
located inside the hollow sphere in the continuum. 

C. Equilibrium magnetization 

We want to determine the thermodynamic equilibrium relation M(H) between the mag- 
netization 

1 iV 
M=^£<m,} = — < m > (2-10) 

and the macroscopic magnetic field H in the thermodynamic limit. Instead of considering 
the dipolar interaction of all particles in an external field in the statistical mechanical prob- 
lem (|2.10|) we take explicitly only interactions between those particles into account whose 
distance is smaller than the sphere radius R s . The other interactions are represented by the 
far-field continuum approximation H s = H + M/3. The magnetization M sp h ere (H s ) result- 
ing from this decomposition of fields is then identified with the equilibrium magnetization 
M(H) of the ferrofluid 



February 1, 2008 9 

M sphere (H + M/3) = M(H) . (2.11) 

Thus after having obtained the approximate expressing for M sphere as a function of H + M/3 
we then obtain from solving ( |2.11| ) for M an approximation for the sought after equilibrium 



relation M(H). The functional dependence of M on H is independent of the probe geometry. 

In the limit of weak dipolar coupling or when R s becomes smaller than the mean distance 
between the particles we find M sv h ere = M sat C (H + M/3) so that the Weiss model is 
recovered as discussed above. 

The magnetization M sp h ere in fl2.11|) depends on two dimensionless parameters that char- 
acterize the thermodynamic state of the ferrofluid. One of these parameters is the volume 
concentration of the particles 

*=^£. (2.12) 
v V 6 y ! 

The ratio <p m ag of the volume of the magnetic material to the total volume is ma9 = 
{Dmagl D) z cj). The other parameter is 

(2.13) 



4:iin kTD 3 ' 

the ratio between a typical energy of dipole-dipole interaction of particles in contact (i. e., 
at the distance of the hard core diameter D) and the thermal energy kT. 



III. CANONICAL PARTITION FUNCTION 

We use the canonical ensemble average to evaluate M sp h ere . Given a system of N inter- 
acting particles with an interaction potential V4, (1 < i, j < N) and external potential per 
particle Vi, the canonical partition function is given by 

Z = J e ~^ Vk ~^<i Vii dY . (3.1) 

Here v% = Vi/kT, Vij = Vij/kT, and dT means integration over the configuration space. In 
our case, a configuration is characterized by specifying the position vector Xj and two angles 
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for each magnetic particle. The two angles define the direction of the magnetic moment 
rrij. The modulus m is assumed to be constant and the same for all particles. Note that 
we ignore any translational and rotational degrees of freedom of the particles that carry the 
magnetic moments, since they have no effect on the magnetization. Only the locations of 
the moments, i. e., of the particles and the orientations of the moments are considered as 
statistical variables. 

In the first -principles statistical mechanical problem identified by a superscript 0, the 
external potential would be the energy of a dipole in the external magnetic field 

V? = -m« • H e . (3.2) 

The interparticle potential is modelled by a dipole-dipole (DD) interaction plus hard core 
(HC) repulsion. Thus 

VS = + V« C , (3.3) 

v o,dd = 3(m^r^)(m j -f^)-m^m i 
ij 4717x04 

for rn > D 

VS C =\ 13 . (3.5) 

00 for r\j < D 

The replacement of the dipolar magnetic fields from far-away particles by a field that 
has its origin in a magnetic continuum results in a new canonical partition function with 
the "external" potential of a dipole 

^ = -m 4 -H s , (3.6) 

in the field H s = H + M/3 and a dipolar interaction term with a cutoff, i. e., 



V 

13 



DD 



V% DD for m < R s 
for r,ij > R s 



February 1, 2008 11 

Note that when using ( |3.6| ) and (|3.7|) in the expression ( |3.1|) for the partition function 
we describe every particle i as being at the center of a sphere of radius R s inside a magnetic 
continuum such that each particle feels the "external" field H s and explicit dipolar fields 
Hdi P oie,near(xi) ( |2 . 7| ) from the particles whose distance is smaller than R s . 

The aforementioned statistical mechanical problems with the long range nature of the 
bare dipolar interactions is thus circumvented by the cutoff at R s in ( |3.7|) that results from 
decomposing dipolar fields into a near and a far contribution. Dipolar forces appear explicitly 
only as forces with a finite range. Their influence on the magnetization in our approach is 
therefore independent of the geometry of the probe. The geometry dependence enters only 
via the effective "external" field H s from the far field contribution. 



A. Born Mayer expansion method 

Since an integral such as ( |3.1|) is hard to solve even numerically we use the Born-Mayer 



expansion method |27j to get analytical results. The key point of this method is to write 

z= jn^na+^r , (3.8) 

k i<j 

where 

fij = e"^ - 1 . (3.9) 

If the typical interaction energy is small compared to kT, the fij can be considered as small 
parameters, and Z can be expanded into a series: 



Z= j Y[e Vm dT + 

m 



i<j k<l 



(3.10) 



These integrals can be factorized and are therefore easier to handle. 
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The first order of ( 3.10|) contains terms like 

'Y[ e - Vm e- V ?* D - V ?° dldtt . (3.11) 



Here d x d Q, is an abbreviation for dxi...dxjvrf^i..<ifijv; and dVti means the integration over 
the possible orientations of m,;. 



B. Expansion in powers of v DD 

Obviously even the first order still cannot be evaluated analytically. Therefore a second 
series expansion is made 

k = e-^V^ - 1 = /<?> +/«+/£> + (3.12a) 



/S° = (e-^-l) (3.12b) 
# } = ^J- e -^°... , a >l . (3.12c) 



ft! 



So we expand /y in powers of the reduced dipolar interaction v® D . The integrals in (|3.10 ) 
that remain to be solved are of the form 

U e ~ Vm fifA? ] - dld ^ ■ ( 3 - 13 ) 

m 

We now introduce a modification of the common graphical representation of Born-Mayer 
integrals as follows 

1. Every distinct particle that appears via interaction terms of the form is represented 
by a circle. 

2. A zeroth order interaction term is represented by an overlap of the circles % and j. 

3. First, second, ... order interaction is represented by one, two, ... lines connecting the 
circles. 
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Note that the representation of the zeroth order dipolar interaction by two overlapping circles 
is a reminder that in this case the integrand is nonzero only if the particles are assumed to 
be in a configuration in which they would indeed overlap. 

It turns out, that the expansion in terms of the f-j means an expansion in powers 
of the two parameters e and 0, that define the thermodynamical system. Every line in a 
representing graph, i. e., every power of v® D results in a factor e. Every n-particle subgraph 
in which all circles are connected to each other directly or indirectly gives a factor of 0™" 1 . 
In the next two sections we will present two expansions considering different terms. 



IV. EXPANSION UP TO FIRST ORDER IN 

In this section only terms up to 0(0) will be taken into account. 



A. Partition function 



In 0(0) the canonical partition function reads 

z = I II e ~ Vk dxdQ + 

J k 

j\{e~"*Y.kd*dn+0(4> 2 ) . (4.1) 

k i<j 

The fij have yet to be expanded in powers of v® D . Figure § shows the corresponding graphs. 

T2 



There are N(N — l)/2 w N /2 ways to choose i and j. Because all particles are identical 
one can write 



Z = f ]Je r w * dx dn + 

jU e ~ Vk f^d^dn+0(<j) 2 ) . (4.2) 



k 

~2~ 



k 

Integrating over most degrees of freedom results in 

N 2 



7V f 

+O(0 2 ) . (4.3) 
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Here 

~ m sinh a s , , , s 

^o = ^o ; z = AnV (4.4) 

a s 

is the partition function of a paramagnetic gas of noninteracting particles in the field H s 
defining the Langevin parameter 

mH s 

a, = V ' < 4 ' B ' 



B. Expansion in the dipolar interaction 



Now we expand fi2 appearing in ( f4.3|) in a power series in e. The n-th summand of this 
series contains integrals of the form 

A n = J e-^ 2 ^ dxjdxadfii^ • (4.6) 
A is special. Here one gets 

Ao = U 8 ^p-Y J (e-5 c - l) dx t dx 2 (4.7) 



The integrand vanishes if r±2 > D. Otherwise its value is —1. Thus 



4 D 3 9 , . 

= "3*1^° ' ( } 



or by expressing the result in terms of 



A = --# 2 . (4.9) 



For n > 1 we have 



n! 

xt-Cj^'^^i^ffii^ . (4.10) 
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Switching from r 2 to the relative coordinate r 12 and integrating over gives a factor of V. 
Then r±2 runs over the sphere volume. We introduce spherical coordinates, i. e., 

mi = m(cos <pi sin ^i, sin (fi sin^i, cos^i) 

m 2 = m(cos if2 sin ^2, sin ip 2 sin $ 2 , cos $2) (4-H) 
r i2 = r i2 (cos ip sin sin sin 1?, cos €) . 

The direction of the magnetic field defines the z-axis. Then, the integral assumes the form 

V r 

— _ g«s COstfl+Qs COS $2 



nl 



2 

m 



xe"" 12 r\ 2 dri 2 du 12 drLidVL2 . 

The new spherical angle CJ12 represents tp and i9. The exact form of the function P is not 
important, but P and therefore P n is a polynomial in the cos and sin of the six angles. 
Integration over four of them can be done analytically. Finally this can also be done for $1 
and $2 by substituting ui )2 = cos$i i2 - One gets an expression of the form 

^l G "^ C (±Sf£) r » dr » ■ (4i3) 

Here we have introduced the correct bounds of the last remaining integral explicitly. By 
setting the lower bound to D we have incorporated the hard core factor. The upper bound 
is given by the cutoff radius R s for the near-field dipolar contribution. While the evaluation 
of G* n can be done analytically, it is quite difficult to do this by hand even for n — 2. We 
therefore used the computer algebra system mathematica to perform the integrations. See 
Appendix |A| for the form of the G* n . 

For n > 2 one can safely set R s = 00 (see below). For n — 1 this would result in a 
logarithmic divergence of the integral. But G\ = anyway, because the calculation of G\ 
involves an averaging over a dipolar field. So by using e and one finally has 



2V 2 
Nn(n- l)n\ 



A n = W7Z7Z 7^ G nMe n cf> n > 2 , (4.14) 
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A l = . (4.15) 

Note that A\ vanishes only in our spherical configuration with finite R s . The divergence of Ai 
in the general, spatially unrestricted case is just an expression of fact that the dipolar forces 
are long range. By treating the distant parts of the ferrofluid as a continuum we incorporate 
any long-range effects and the resulting geometry dependence via the field H s = H + M/3. 
Into this field enters the relation between the external and the macroscopic internal field. 

Two further comments should be made here. A generalization of our calculation for 
central symmetric interactions other than a hard sphere potential is possible. It requires 
an analytical or numerical evaluation of integrals of the form / r 2 ~ 3n e~ v R dr in (|4.13|) , with 
v sr _ ySRj fcji an j ySR denoting the r-dependent short range potential in question. 

The second thing is that we can now make quantitative statements about how large 
the virtual sphere has to be chosen. To ensure A\ to vanish unambigously in (|4.13|) and 
to introduce H s instead of H e as "external" field R s has only to be finite. The larger R s 
the better is the modeling of large-distance particle correlations entering into A n for n > 1. 
Taking the limit R s — > oo as the final step in the calculation of the A n is therefore appropriate 
from this point of view. On the other hand, the requirement of uniformity of the fields H 
and M, that allows us to write H s = H + M/3 restricts R s to values below the scale on 
which H and M vary. If one would use a finite radius R s one would get instead of (|4.14j ) for 
n > 2 

2V 2 



N7r{n — \)n\ 



(D/R s f n - S (4.16) 



which allows an error estimate: Consider a system where M and H do not vary on the scale 
of, say, /im. For ferrofluids, D ~ 10 nm. Choosing R s = 10D or 100D is then both allowed 
and implies a difference in A 2 of about 0.1 percent. The result for R s = 100D is better than 
for R s = 10D, because in the latter case particles in a distance range between 100 nm and 
1 fim are treated in the continuum approximation and not correctly. But the error that is 
made by treating the ferrofluid as a continuum already beyond R s = 10D is only about 0.1 
percent. We can safely assume that the macroscopic, magnetic properties do not vary on 
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this scale. Thus 100 nm is an appropriate medium scale on which both requirements hold: 
The continuum approximation works well beyond this cutoff radius and the macroscopic 
fields H and M should be constant on this scale. Except for the calculation of Ai, it is then 
possible to set R s = oo in the calculations of the integrals. 

Using the results ( |4.14| ), ( |4.15|) , (|4.9|), and ( [4.6| ) in ([4.3[), one gets the following expression 
for Z: 



Z = Zn 



l-4N(j) + N<f)^rG n (a s )e r 



n=2 



+ O(0 2 ) . 



Here we introduced the functions 



G n (a s ) 



( a s 



167r 3 (n — l)n! Vsinha; 
some of which are given in Appendix |A]. 



G*n( a s) , 



(4.17) 



(4.18) 



C. Free energy and magnetization 



The next step is to compute the free energy 
F 



kT 



-InZ = -iVlnzo 



-In 

+0(0 



l-AN<j ) + N ( f)J2G n (a s )e n 



n=2 



(4.19) 



In 0((p), we can use ln(l + x) = 1 + x here: 



kT 



-N In z + 4iV0 -N<f>J2 G n (a s )e r 



n=2 



+O(0 2 



(4.20) 



The magnetization turns out to be 



M, 



1 OF 



m OF 



sphere 



a. 



HoV dH s [i VkT da s 



fM V 

+ou 2 



n=2 



(4.21) 
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The leading term is the Langevin function £ times the saturation magnetization M sat = 
of the fluid. 



In order to determine M(H) we identify, according to ( |2.11|) , M sphere (ot s ) with M(a), 
i. e., 

M „ / mM\ 



£ a + 



M sat V 3kT, 

+0E^(« + ^) e " + O(0 2 ) , (4.22) 

where a is the usual Langevin parameter, 

mH 

a = w ■ (423) 

Instead of trying to find the function M(a) that solves this equation exactly we expand the 
functions £ and G' n for small into a series around M = and reinsert it on the right hand 
side. Using the fact that m ^ at = 8(pe grows linearly in <p we arrive at 

= L ,o + £ L liB c" + 0(4> 2 ) (4.24a) 

M sat n= l 



with 

L 0i0 = £(a) (4.24b) 

Li,i = 8£(a)£'(a) (4.24c) 

L l n = G' n (a) for n > 2 . (4.24d) 

This is a consistent approximation in terms of 0. On the other hand, solving Eq. (|4.22|) in a 
formally exact manner for M would introduce higher orders in <p which we already neglected 
to arrive at Eq. (fL2^) . 



Note also that M sp h ere (a s )/M sat ( |4.21| ) contains explicitly a term ~ 0e 2 as lowest non- 



trivial power coming from the expansion in the near-field dipolar coupling strength. On the 
other hand the self-consistent solution ( |4.24| ) that solves Eq. ( f4.22| ) starts out with a contri- 



bution ~ <pe. The latter arises from the far-field dipolar continuum via the magnetization M 
in the dipole- induced shift of the argument, a + T^f, of the Langevin function in Eq. (|4.22|) 
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- in the absence of any dipolar interactions in the system one would have H s 
leading to ideal paramagnetism. 
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= H e = H 



D. Comparison with previous results 

The Onsager model, the Weiss model, and our calculation agree that up to order 0e 

U = £(a) + 80e£(a)£'(a) + 0(</> 2 ) + 0(e 2 ) . (4.25) 



M, 



sat 



This expression was also derived by Buyevich and Ivanov |22] with a calculation similar to 
ours. However, they did not introduce a magnetic continuum approximation. Instead, they 
assumed a special probe geometry of a long cylinder parallel to the external magnetic field 
and performed an integration over all the particle's dipolar fields in the cylinder explicitly. 
The magnetization was therefore given in terms of the external field. Their result agrees 
with ours because for the cylindrical geometry chosen in |22] H e equals H. 

A second paper that deals with our problem in a similar way was published by Kalik- 
manov [[53J]. In section 4, the author arrives at an equation for the magnetization that reads 
in our notation 

= £(«) + 3^G' 2 (a) r 9 -^dx . (4.26) 

M sat J I X 

Here go(x) is the hard sphere correlation function. In our (^-approximation this function 
has to be set to one. Then the 0e 2 -term agrees with ours. Note, however, that the above 
result (|4.26 ) of Kalikmanov does not contain the 0e-term resulting from the magnetic field 



from the continuum. 



V. EXPANSION UP TO SECOND ORDER IN <p AND e 

It is possible to calculate 0(</> 2 )-terms of the Born-Mayer expansion when e is taken into 
account up to second order only. A more elegant way to calculate the magnetization in this 
approximation makes use of the grand canonical rather than the canonical ensemble. This 
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approach allows to avoid the determination of some terms that can be factorized into already 
known integrals and cancel out in the calculation of the free energy. However, the grand 
canonical approach has the disadvantage that it yields the magnetization as function of the 
chemical potential fi rather than the particle number N. Some more algebra is then required 
to find out the function fi(N). Here we continue to work with the canonical ensemble. 



Figure ||] shows the 12 additional graphs that are of second order in and of less than 
third order in e. Four of them vanish because they contain at least one first-order dipolar 
interaction term between otherwise unrelated particles. Integration over the relative position 
of these particles while leaving the relative positions between all other particles and the 
direction of the magnetic moments fixed yields zero since it involves a spatial averaging 
over a dipolar field. The graph labelled with the letter F vanishes for similar reasons that 
are explained in Appendix [FJ where we calculate the integrals one by one. Their respective 
contribution to the partition function is 



A. The graphs 




32A0 2 



(5.1a) 




16N(f) 2 e 2 G 2 {a s ) 



(5.1b) 




8(A 2 - 6N)(f) 2 



(5.1c) 




4(A 2 - 6A)0 2 e 2 G 2 (a s ) 



(5.1d) 




5N(j) 2 



(5.1e) 








(5.1f) 



l + 61n2 

4 



N(j> 2 e 2 G 2 (a s ) 



(5.1g) 




N<j) 2 e 2 K(a s ) 



(5.1h) 



The functions G 2 and K are given in Appendix [A]. 
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B. Free energy and magnetization 

Now we have all necessary terms at hand to calculate the canonical partition function 
up to the desired order: 

^-=l-4(N-l)(f>+(N- l)0e 2 G 2 (a s ) 

+32iV0 2 - 16N(j) 2 e 2 G 2 (a s ) + 8(N 2 - 6iV)0 2 

-4(iV 2 - 6AO0 2 e 2 G 2 (a,) - 5iV0 2 (5.2) 
+ l + 6 A ln2 N ( p 2 e 2 G 2 (a s ) - N<p 2 e 2 K(a s ) + H.O.T. . 



The terms in 0(<f)) appear already in (|4.17 ). They are presented here including the next 
higher order in N. The other terms come from - Z^. To include all terms of 0(4> 2 , e 2 ) in 
the free energy one has to approximate the logarithm ln(l + x) by x — x 2 /2. The quadratic 
order is necessary only for the 0(<f>)— terms. New terms of 0(N 2 ) appear and cancel against 
those of the terms from Zc and Zr>. One gets 

— = -N In z + AN<p + 5A^0 2 
kT 

-Nct)e 2 G 2 {a s ) - 1 + 6 ^ N<p 2 e 2 G 2 {a s ) (5.3) 
+N(j) 2 e 2 K(a s ) + H.O.T. . 

The result is proportional to N as it has to be. 
The magnetization of the sphere is 

Msphere(as ± = £(a s ) + <Pe 2 G' 2 (a s ) (5.4) 



M sat 

+ 1 + 6ln2 <P 2 e 2 G 2 ( as ) - <P 2 e 2 K'(a s ) + H.O.T. 



To calculate the magnetization as a function of a we identify (|5.4j) with M and use again 



a s = a + The right hand side of ( |5.4j ) has now to be expanded around a up to second 



order and the resulting equation has to be iterated twice to take into account all important 
terms up to e 2 2 . The result is 
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M(n\ 

= L ,o + <j>eL lt i + 4>e L lj2 + (ft e L 2j2 + ••■ (5.5a) 



M sat 

with Lo,o,-^i,i, and Li >2 denned in eqs. (|4.24b| ) - ( |4.24d| ) and 



L 2 , 2 = 64£(a)£'(a) + 32£(a)£"(a) 

1 1 + 61n2 -G' 2 (a) -K\a) . (5.5b) 



4 

For the discussion in the next Sec. we decompose 

L 2j2 (a) = L s 2 f ere (a) + L 2 te 2 rative (a) . (5.6a) 

The function 



7- sphere 1 -|- 6 In 2 / . 

^2,2 — 1 ^2 — A 



(5.6b) 



occurs already in the expression ( |5.4| ) for the magnetization M sp h ere (a s ) of the sphere. The 
contribution 

L iterative = 64 £(£')2 + ^CC" (5.6c) 

arises in obtaining the selfconsistent solution of the equation M = M sv h ere with an expansion 
and iteration. 

VI. DISCUSSION OF THE RESULTS 

We will first show that our result ( |5.4|) for M sphere (H s ) does not lead to a ferromagnetic 
solution in contradistinction to the Weiss model. Then we discuss the behavior of the 
different terms contributing to ( |4.24| ) and to ( |5.5| ) and we delineate the range of reliability 
of the simplest approximation. Finally, we address problems arising when comparing with 
experiments. 
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A. Spontaneous magnetization? 

Investigations based on density functional methods by Groh and Dietrich |4| and on 
Monte Carlo methods by Weis and Levesque |8||J provided support for the existence of 
magnetized phases for absent external field H e , i. e., ferromagnetism, in the system of 
dipolar hard spheres we consider in this work. Groh and Dietrich consider a ferrofluid 
probe of needle-like shape where H = H e and find a transition to a magnetized phase at 
<pe « 0.35. But they consider this value as being overestimated and refer to ||. Weis 
and Levesque study a case without demagnetizing fields, i. e., again H = H e . They find a 
transition to a magnetized phase at e = 6.25 for <fi ps 0.35. As discussed in detail below, 
these values are outside the range of reliability of our results. 

The Weiss model does also show ferromagnetic behavior. It is recovered from (|5.4j) by 
keeping only the leading-order term C(a s ) describing a single moment in the field H s = 
H + M/3. The resulting self-consistency equation 

allows for zero field a solution with finite magnetization when kT < mM sat /9. Using fl2.4| ) 
combined with fl2.12p and ( |2.13j ) this condition is equivalent to <pe > 3/8, about the same 



M = M^Zt + y ) = M sat C 



value as in ||. So according to the Weiss model the ferrofluid will show ferromagnetic 
behavior below a critical temperature that grows linearly with the saturation magnetization 
M sat of the ferrofluid. But even for a ferrofluid consisting of cobalt particles with a magnetic 
core diameter of 10 nm and a magnetic volume fraction of <p m ag — 0.1 the critical temperature 
would be as low as 90 K. 

While the transition combination e = 6.25, <fi ~ 0.35 of is outside the range of 
reliability of our results, the threshold location <pe = 8/3 of the Weiss model may be not. 
However, in agreement with we do not find selfconsistent ferromagnetic solutions of the 
equation ( |5.4j ) M = M sp h ere (H + M/3) within this range. We have numerically confirmed 
that for H = the equation M = M sp / iere (M/3) allows always only the trivial solution 
M = 0. 
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B. Contribution from different orders 



Now we will take a closer look on the functions of a involved in ( 4.24 ) and (5.5). All these 
functions are odd as it has to be for reasons of symmetry. For a — > oo they vanish as I /a 2 
or faster. Because 1 — C(a) ~ 1/at that means that the predicted magnetization is always 
smaller than the saturation magnetization for a — > oo. Nevertheless the magnetization 
can assume unphysical values > M sat for intermediate a if e or (ft is big enough for the 
approximations to become invalid. 



1. Behavior in linear order of < 



We will first discuss the result ( [4.24|) for the magnetization that was obtained up to 
linear order in the volume fraction cf). In figure |] the functions L\ % \ and L 12 are plotted. 
The values of the higher-order functions are smaller, but their shape remains more or less 
the same as the logarithmic plot in figure [| shows. Because L\^ n and Li, n +2 differ by about 
one order of magnitude one can conclude that by including higher and higher orders of e the 



series (5J3) for the magnetization converges, as long as e is smaller than pa 3. For this large 
value of e strong agglomeration can already be expected. 

For small a, L l n is proportional to a (a 3 ) for odd (even) n. The initial susceptibility 
can therefore be written as 



X(H = 0) = 
Xo(H = 0) 



l + 0£ s Wie 2n+1 + O(0 2 

n=0 



(6.2) 



Here 



Xo(H = 0) 



mM, 



sat 



3kT 



(6.3) 



is the initial susceptibility of the ideal paramagnetic gas, and the nonvanishing si ra we 
calculated are 



February 1, 2008 25 

8 8 32 



3 ; *i,3 ?5 ; s li5 3g75 
8 148 
Sl ' 7 ~ 19845' S1,9 ~ 12006225 ' ^ ' ' 



Figure |^ shows Xo(H = 0) (thick dashed line), and the susceptibility x(H = 0) 
including progressive orders <f>e, 0e 3 , 0e 5 , 0e 7 , and 0e 9 (thin dashed lines, from bottom to 
top) as a function of e for <ft = 0.15. The sequence of these thin dashed lines shows that 
this series converges in the e-range of figure The last thick full line in Figure |] represents 
x{H = 0) including the contributions in order 2 e 2 . It shows that the latter are even for 
= 0.15 not yet important. 



2. Behavior in second order of 4> 

Now we take a look at the functions L s ^ ere (|5~6|b) and L^{ aUve Q5"^c) that add up to 
-^2,2 ( p.6| a) which enters in order <fr 2 e 2 into the magnetization ( |5.5| a). 

Figure |^ shows that the contributions L S 2^ ere and Jj^ atme almost cancel each other at 
small a. This is why the influence of the 2 e 2 -terms on the susceptibility in figure |6| is so 
small. However, at higher a the 2 e 2 L 2) 2^term becomes important. Comparing the latter 
with the linear one, (peLi^, one finds that they contribute for e<p ~ 0.5 equally at larger a. 

Except for very small a £2,2 is negative, because it includes higher-order particle position 
correlations that result in a better modeling of the distance distribution due to the finite 
size of the particles. The mean distance is bigger in this approximation and the induced 
dipolar fields at the particle positions are therefore smaller. 

The influence of the 2 e 2 L 2j 2 contribution to the magnetization is shown in figure |8] for 
e = 2 and <p — 0.05. For these parameters this term is already large enough to cancel 
almost exactly the sum of all contributions Li ra 0e n , with n > 2 from the linear order in <p at 
moderate a. Figure [5] shows the susceptibility x{H) = dA g^ f° r the same parameters. At 
higher a, the cancellation of the higher L l n -terms against the /^^-contribution can again 
be seen. At smaller a, however, the behavior is different. There the contribution of the 
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L ln -terms is much larger, whereas the L 2j 2~contributions vanish. 
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C. Reliability of the 0(4>e) approximation 



We can determine the range of reliability of the simplest approximation 

M 



M, 



sat 



L , (a) + Li,i(a)0e 
£(a) + 8£(a)£'(a)0e 



(6.5) 



to the magnetization that includes effects of dipolar interactions since we know the higher- 
order corrections in as well as in e. To that end we investigated the ratios 



O(0e n )-terms (n > 1) 



and 



L Qfi (a) + £i,i(a)0e 



O(0 2 )-terms 



(6.6) 



(6.7) 



L 0)0 (a) + £i i i(a)0e 

The first ratio assumes its maximum at a = 0, that means the initial susceptibility is most 
sensitive to higher order corrections in e. The second ratio (|6.7|) assumes its maximum 
around a = 2, that is near the maximum of the absolute value of the numerator (as seen in 
figure ^). 

In the e~4> plane of figure 0(a) we show isolines of the maximal - with respect to a - 
ratio (6^) and figure 0(b) shows the analogous isolines for the ratio ( |6.7|) . The comparison 
shows that the smallness of e is more important to keep the ratio (|6.6j ) small, whereas in 
( |6.7|) the value of (f) is also important. As rules of thumb one can say that the approximation 
( |6.5| ) is valid within about 1-2 percent if e < 1 and e<p < 0.04. If the first constraint is not 
fulfilled, higher orders in e have to be taken into account. Higher orders in <p are needed if 
the second constraint is not fulfilled. 
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D. Comparison with experiments? 



There are several papers 23-|2q| that aim at investigating the influence of dipolar in- 



teractions on the magnetization by comparing theoretical models developed so far with 



experimental magnetizations of ferrofluids. The mean spherical model [21] was reported 



to show good agreement with experiments. Pshenichnikov [23|] found also good agreement 



with the high temperature approximation [P2"| , i. e., the approximation (|6.5|). But this ansatz 
failed in the magneto-granulo metric analysis done in ||26|| . 

We do not present a comparison of our results with the experiments on the magnetization 
in the literature because of several problems: In our theory it is necessary to distinguish 
between the particle diameter D and the magnetic core diameter D mag that is found in 
magneto-granulometric measurements. This problem does not arise in the mean spherical 
model or the high temperature approximation, where e and <fi enter only via the factor 
( t )e = 2W™okT = M 2AkT ^ na t * s independent of D. Also, corrections such as the temperature 
dependence of the saturation magnetization or the fluid density should be taken into account 



But the major problem in comparing directly with experiments is that our theory does 
not take into account the polydispersity of ferrofluids. The effect of polydispersity is signifi- 
cant already in the absence of any dipolar interaction. This can be inferred from the dashed 
and the dotted curves in Fig. |ll] representing the reduced magnetization of noninteracting 
magnetic particles having a polydisperse and a monodisperse distribution of particle diam- 
eters, respectively. Here the common particle diameters of the latter is .D 3 ^ 3 , where D 3 is 
the third moment of the particle size distribution 

1 In 2 D/D a 

of the former. The mean magnetic moment m and the saturation magnetization of the two 
systems are the same. For comparison with the effect of dipolar interaction in monodisperse 
systems the full curve in Fig. [H] shows our result for M ( |5.5| a) including all terms ~ <pe n 
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and the term <p 2 e 2 . Hence the effects of polydispersiveness alone, i. e., without interaction 
are comparable in size with the effect of dipolar interactions in monodisperse systems. Thus 
clearly an extension of the here presented Born-Mayer expansion method to the case of 
polydisperse interacting particles is desirable. 



VII. CONCLUSION 

We calculated the free energy and in particular the magnetization M of a ferrofluid 
as a function of the macroscopic magnetic field H . To do so, we used the technique of 
the Born-Mayer expansion together with an expansion in terms of the dipolar coupling 
energy. The magnetic particles were assumed to be hard spheres with a common hard core 
diameter D and magnetic moment m that interact via long range dipolar interactions. This 
feature may result in a geometry dependence of thermodynamic properties. We treated this 
problem by dividing the dipolar field at some position that is produced by the magnetic 
moments of the particles into a near-field and into a far-field part depending on whether the 
particle distance from Xj is larger than some radius R s or not. In this way every magnetic 
particle is imagined to be located in the center of a sphere of radius R s . The far-field 
dipolar contribution from particles beyond R s is then replaced by a magnetic continuum 
with magnetization M and magnetic field H. Here R s is chosen to be such that M and H 
are homogeneous on the scale of R s . The magnetic continuum outside the sphere produces 
in the center of the sphere the magnetic field H s = H + M/3. This field acts as an "external" 
field on the particle in the center of the sphere. The near-field interaction of the latter with 
the other particles within the sphere being at a distance smaller than R s is treated explicitly. 
Thus in our statistical mechanical calculations there appear dipolar interactions only with 
interparticle distances less than R s . However, since the cutoff dependence of the relevant 
expressions occurring in these calculations is negligible already beyond a radius of the order 
of 10D ~ 100 nm we used R s = oo in these expressions. 

The expansion of the partition function for these interacting particles in terms of the 
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volume ratio and the dipolar coupling strength e yields an expression for the magnetization 



as a function of the "external" part of the field inside the sphere. The magnetization M sp h ere 
is then identified with the magnetization M(H) inside the continuum so that a selfconsistent 
relation results. The aforementioned geometry dependence of M in the general case is 
incorporated via H. 

We presented two different expansions in e and 0, one containing only linear terms in 
<f), the other also second order <fr terms, but only up to 0(e 2 ). We discussed the range of 
applicability in the 0-e plane of their results for M(H) and compared them to the most 
simple approximation to the magnetization that contains the dipolar effects only in linear 
order in e and 0. The selfconsistent relation for M(H) that contains only up to second 
order terms in both parameters does not admit a ferromagnetic solution with spontaneous 
magnetization. Finally we showed that an extension to polydisperse interacting particles is 
desirable. 



M sp h, 



ere 



M sphere (H + M/3) 



(7.1) 



ACKNOWLEDGMENTS 



This work was supported by the Deutsche Forschungsgemeinschaft (SFB 277). 



APPENDIX A: THE FUNCTIONS G n AND K 



The functions G* n (x) in ( f4.13|) are related to G n (x) via ( }4.18|) : 




(Al) 



The functions G n {x) introduced in ( |4.18| ) have the form 




(A2) 
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where the functions G$(y) are polynomials. The first four triple are given by 

G? ) (y) = | + |y a + yy 4 (A3*) 

G? ) (y) = -|y-yy s (A3b) 

Gf (v) = y l/ 2 (A3c) 



cS 0) (y) 
^(y) 

Gi a) (y) 



£ 2 _ 48 4 _ 12 6 
35 y 35 y 7 27 



35 
16 

105 



V 



8 3 24 - 

-y H u 

5 y 7 y 

8 2 12 4 
35* ~T V 



(A4a) 
(A4b) 
(A4c) 



™ = W 5 + 35^ + 35^ + y^ 6 + 122/8 ^ 
Gi 1} (y) = -^y - jjy 3 - yy 5 - 24y 7 (A5b) 

^(v) = ^y 2 + ^WlV (A5c) 



, l/ x 12 2 208 4 852 fi 480 8 540 in /A , 

G^Hy) = y 2 y 4 y 6 y 8 y w A6a 

5 vy; 385 y 385 y 77 y 11 y 11 y V ; 

-m, , 8 16 3 472 5 600 7 1080 9 , . 

G = fe) = + 385* + 77* + ~n y + — » (A6b) 

G f\ y) = + JL? _ « * _ H y « _ 5iV • (A6c) 
5 yyi 1155 1155 77 11 11 V ; 

All functions G W (x) have a well defined limit for x — > although this is not obvious for the 
above explicit expressions. Their values at a; = are closely related to the coefficients in 
the e-expansion of the second virial coefficient for the system of dipolar hard spheres in the 
absence of a magnetic field. The calculation of this coefficient dates back to |H| and can 
also be found in p6| . 

The function K ( |B22| ) that appears in the O(0 2 )-terms of the free energy is given by 
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K(x) = — coth 3 x + ( — + 12 ) coth z x 
X \x 2 



18 24 \ 6 12 /4r7 , 

— + — cothx + — + — . (A7) 



X 3 X J X A X 2 



APPENDIX B: GRAPHS IN SECOND ORDER OF 

Here we determine the contribution to the canonical partition function from the graphs 
A-H shown in Fig. [| There often appear hard core interaction terms that are just expres- 
sions of the requirement that two particles have to or must not overlap. We define two 
abbreviations: 

e -r _ i = _ 0ii , (Bl) 
e'< C = O l3 . (B2) 



1. Graph A 



The graph A represents fijfj^- There are N 3 ways to choose the constituting particles. 
But because j and k are equivalent only A 3 /2 distinctive graphs remain. Integration over 
all variables except the positions of particle j and k relative to i yields 



A3 

IT 

N 3 ( sinha 



2 



i 



The remaining integral factorizes and we can make use of the results for Aq (|4.9|) . The 
contribution of graph A to the partition function is 

Z A = 32NZ (p 2 (B3) 

where Z Q is given by (|4.4|). 
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2. Graph B 

The graph B represents /J^/Jk ■ All three particles appear in different ways, thus there 
are N 3 different graphs. After integration over the degrees of freedom of the noninvolved 
particles and switching to relative coordinates with respect to particle % the integral factorizes 
again and one can make use of the results for A ( |4.9|) and A 2 ( 4.14 ). We get 



Z B = -lMZ 0( j)'e 2 G 2 (a s ) . (B4) 



3. Graph C 

The graph C represents fjj fj^ . Here we have also to include the next higher order term 
when we calculate the number of combinations to get the 0(A r )-terms in the final result: 
There are (N 4 — 6iV 3 )/8 different terms. The integral for graph C can be factorized so that 

Z c = 8(N 2 - 6N)Z (f) 2 . (B5) 

4. Graph D 

The graph D represents fij fj^ ■ The calculation is similar to the calculation of graph 
C. Again, we need the next higher order term in N. There are (iV 4 — 6iV 3 )/4 combinations, 
twice as much as for graph C because the pairs and (k, I) are not identical. One gets 

Z D = -4(iV 2 - 6N)Z (f> 2 e 2 G 2 {a s ) . (B6) 



5. Graph E 

The integral containing the term fifffkfki ls the first really new integral. It involves 
only hard core interactions and does not contribute to the final expression for the magneti- 
zation. But for completeness we will calculate it also. The trivial integrations yield 
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N 



Z E = ~ Un^^ ) V N - 2 I 12 13 23 dr 12 dr 13 



(B7) 



We keep the distance r 12 fixed. The center of particle 3 has then to be inside two spheres 
of radius D around particle 1 and 2. Integrating over the position of particle 3 yields the 
overlap volume V of the two spheres 



V„ = -nD* 



1 _ ^ r j2 J_ ( r u\ 

AD 16 V DJ 



(B8) 



Therefore 



iV 3 Z t 
Ze = — —-^ J Oi2V {ri 2 )dYi2 

N 3 Z rD 



(B9) 



6 V 2 

Performing the last integration results in 



f 

4vr / V (r 12 )rl 2 dr 12 
Jo 



Z R = -5NZ n 



(BIO) 



6. Graph F 

The graph F represents flff^ffk- As already stated this integral vanishes which can be 
seen as follows: Consider an arbitrary configuration belonging to some value of the integrand 



e 



vi-vj-vk AO) AO) Ai) 



WWW ■ (BH) 



While leaving the direction of the magnetic moments fixed the whole configuration can be 
freely rotated around particle j changing only the /j^ term. Integration over the resulting 
configurations involves again an averaging over a dipolar field on a spherical surface. 

7. Graph G 

The calculation of the N 3 /2 integrals belonging to fif'Wffj) is similar to the calculation 
for graph E. First we integrate over all degrees of freedom except the distance between 
particle j = 1 and k = 2 and the position of particle i = 3: 
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x J 13 23 12 r^dr l2 dv 3 . (B12) 
Integrating over r 3 results again in an overlap volume term: 

AJ~3 y, f _ 

Z G = -Y^G 2 (a s )e 2 D e j O l2 V {r l2 )r^dr l2 . (B13) 

Here the lower integration boundary is r 12 = D because of the remaining hard core factor. 
The upper integration boundary is r i2 = 2D because the possibility that particle 3 overlaps 
with particle 1 and 2 is still required. The final result is 

Z G = 1 + 6 A ln2 NZ <j) 2 e 2 G 2 (a s ) . (B14) 



8. Graph H 

The last graph H is the most complicated one. It represents the term fjj^fikfjk that 
appears N 3 /2 times. The problem here is to fulfill the requirement that particles j and k 
have to overlap in terms of properly chosen integration limits. We start with performing the 
trivial integrations: 

Z B = -\n Z z»-"V 

x j e-^-^v° D v° D 23 12 13 (B15) 
xdr 12 dr 13 d£lid£l 2 d£l 3 . 

Whether the integrand vanishes due to the hard core factors depends only on the distances 
r i2, r i3, and the angle $23 between r 12 and r 13 . Consider a special orientation where 

?; 2 = (1,0,0) (Bie) 

f? 3 = (costf 2 ,0,sintf 23 ) , (B17) 
with < -# 2 3 < 7r. A general configuration of the particles' locations can be written as 
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fl(2,3) = 7^(^(tf)7^(V0??(2,3) > (B18) 

where TZ X , lZ y , and 7?. 2 are Eulerian rotation matrices for the angles ip, and (p. Using this 
form the integration over the factors that depend on these angles: 



/•2ir r2n 

/ / / v ?2 D Vn D cos tidipdMip , (B19) 
J o J-tt/2 Jo 



can easily be performed with mathematica. We call the result I(r\2, ri3, $23) nij). 
Next, we integrate over the orientations of the 

J e- Vl ~ V2 - V3 1 (r 12 ,r 13 , ^23, ^d^d^dQs . (B20) 

The result depends only on r 12 , r 13 , and $ 2 3- Using it in ( [B15|) yields 

X /0230 1 20 1 3 m47r2(2COS ^ 23_Sin ^ 23) 



10(47r/x A;r) 2 r 12 r 1 3 
x sin $23drndr 13 d$23 , (B21) 



where 



8 Vsinha s / J-U-i J-i 
x (w 2 + 3)u 2 u 3 duidu2dus ■ (B22) 

The explicit expression for K(a s ) is given in Appendix [A] (eqn. |A7| ). 

Now we discuss the hard core terms. 7*12 and ri3 have to be greater than D to avoid the 
overlap with particle 1. Furthermore \r\2 — < D has to be fulfilled for particle 2 and 3 
to overlap. As a last requirement, $23 has to be smaller than some angle i?^ 8 * that depends 
on r\2 and 7*13. Trigonometry shows that 

r 2 1 J2 _ f)2 

cos ^a X = r 12 + r 13 
2r 12 r 13 

In this configuration the distance between particle 2 and 3 is exactly D. 

We perform the integration over $23 from to f?^ 3 * i n ( P^lp , choose the correct limits 
for r±2 and ri3, and drop all hard core terms: 
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roo pri3+D 

X 



X 



4 2 
WIT 



D Jmm(D,rn-D) 10(47T/! A;T) 2 r 1 2r 1 3 



(B24) 



12 



r 2 3 - D 2 



D 



2\ 3 ' 



2ri 2 ri 3 

The result of the last integration is: 



2ri 2 ri 3 



dr-todr 



12"' 13 , 



3 ^0 



4 2 



— iV : 

3 \/ 2 " v " iV 48(47r/i fcT) 2 
-iVZo0 2 e 2 ^(a s ) . 



(B25) 
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FIGURES 

FIG. 1. Geometry of our model. Every particle (black) feels the magnetic near field generated 
by the dipoles of the neighbours (dark gray) within the radius R s and a contribution from the 
continuum (light gray) that models the fields of the far-away particles and the external field. 



FIG. 2. The first four graphs needed for the expansion of Z in section IV. They correspond to 
the terms f-j (n = 0, 1, 2, 3), and are of the order 4>e n . 

FIG. 3. The twelve additional graphs needed for an 0{4> 2 , e 2 ) expansion of Z . The integrals for 
the crossed out graphs vanish. Graph F vanishes also; see Appendix |B|6. 

FIG. 4. The functions L^i and versus a. Note the different scaling. 

FIG. 5. The functions L\ iU versus a. 

FIG. 6. Initial magnetic susceptibility for <p = 0.15 as a function of e. 



FIG. 7. The weight of the 0(4> )-terms that appear in (|5.5| ) are shown versus a. The terms 
^sphere ^y^) and Jjterative (|0I C ) th & t ao -d U P to L2,2 (|5.6|a) are discussed in the text. For com- 
parison, the O(0e)-term L\ t \ is plotted as well. 

FIG. 8. The reduced magnetization for e = 2 and <p = 0.05 for moderate a. Taking into account 
</>e n -terms results in a higher magnetization than given by the Langevin function. However, all 
contributions from the terms (j)€ n with n > 2 are almost exactly canceled by the contribution from 
the second order term cj) 2 e 2 for the parameters e, <f> considered here. 

FIG. 9. The reduced susceptibility for e = 2 and <fi = 0.05 as a function of a. The higher order 
corrections are largest at a = 0. At moderate a, the cancellation of the terms of order <pe n with 
n > 2 against the term of order (j) 2 e 2 can again be seen. 
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FIG. 10. Quality of the lowest order expression ( |6.5[) for the magnetization, (a) shows the 
isolines of the maximal - with respect to a - ratio fl6.6|) in steps of 0.01 and (b) shows those of the 
ratio flajl ). 

FIG. 11. Comparison of the effects of polydispersity and of dipolar interaction. Plotted is 
the reduced magnetization versus a for different ferrofluid models: a noninteracting monodisperse 
system (only -Lo,o)> a noninteracting polydisperse system, and a monodisperse system with dipolar 
interaction for <j) = 0.05 and e = 2. The polydisperse system has a lognormal distribution of 
particle diameters (eq. |6.8| ) with a typical width of a = 0.3 and the same third moment D 3 as the 
monodisperse fluid. 
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